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ABSTRACT 


A simple method of obtaining rational forms as approximations 
to functions which may be expressed as power series is developed. The 
method is near-optimal under the Tchebycheff norm. While most approaches 
are iterative in nature, the method presented here is free of this 
characteristic. 

The technique is developed as a combination of Padd rational 
approximation and Lanczos'Tau method, which uses Tchebycheff polynomial 
properties for improving accuracy. In addition to some simple non- 
linear functions, some examples from the field of astrodynamics are 
used to illustrate the method. 



ACKNOWLEDGMENTS 


I extend my thanks to Dr. Ray Nachlinger for serving as chairman 
of niy thesis committee and to committee members. Dr. Bart Childs, Dr. 
George -Born, and Dr. Pat Hedgecoxe. 

In particular I thank my wife, Jean, for her encouragement and 
for typing this thesis. 


v 



TABLE OF CONTENTS 


Page 

ACKNOWLEDGMENTS i V 

LIST OF FIGURES viii 

Chapter 

1. INTRODUCTION 1 

1.1 Rational Forms as Approximating- 

Functions 1 

1.2 Background Study of Rational Approximation .... 3 

1.3 Scope of the Investigation 4 

2. APPROXIMATION CHARACTERISTICS 5 

2.1 Existence of Best Rational Approximations . . . . 5 

2.2 Characterization of Best Rational 

Approximations 7 

3. PADE APPROXIMATION 9 

3.1 Development of the Pad£ Method ........... 9 

3.2 Example - Tan’* 1 (x) 12 

3.3 Error Expression 16 

. 4. THE TAU METHOD 18 

4.1 Tchebycheff Polynomials 18 

4.2 The Tau Method 19 

vi 



vii 

5. TAU AND FADE METHODS COMBINED 24 

5.1 Development of the Tau-Padd Equations 24 

5.2 Examples 29 

5.3 Error Expression 37 

5.4 Applications to Kepler’s Equation 40 

6. EXPLICIT FORMS AND GENERALIZATIONS • 55 

6.1 Explicit Expressions 55 

6.2 Generalization 58- 

7. SUMMARY . 62 

BIBLIOGRAPHY 65 

APPENDIX 67 



LIST OF FIGURES 


Figure Page 

2.1 Condition of Optimality 'j - # 6 

3.1 Error Curves, for tan” 1 (x) (Padd Method) 15 

4.1 Error Curves for e x (Tau Method) 22 

5.1 Matrix Form of Tau-Padd Coefficient 

Equations 27 

5.2 Error Curves for e x 32 ■ 

5.3 Error Curves for An(l+x) 34 

5.4 Error Curves for tan _1 (x) * 38 

5.5a Error for Sin x . . . 43 

5.5b Error for Sinh x 43 

5.6a Error Curves for E 44 

5.6b Error Curves for H 44 

5.7 Error for Quadratic Approximation to 

Sin and Cos 47 

5.8 Eccentric Anomaly Error 49 

5.9a Error Curves for C 53 

5.9b Error Curves for S 54 

viii 



CHAPTER 1 


INTRODUCTION 


The concept of rational approximation is introduced and set in 
perspective with the more familiar polynomial approximation. The 
historical development is traced briefly and the scope of the investi- 
gation is set forth. 


1.1 Rational Forms as Approximating Functions 

A common method of approximating transcendental and other non- 
linear functions is through the use of polynomial approximation. It is 
an attractive approach chiefly because of its simplicity, but generally 

requires polynomials of high degree for high accuracy. In recent years 
* 

rational functions as approximating forms have come under investigation. 
Rational functions have been found to offer considerable flexibility and 
accuracy in the approximation of certain functions, the rational form 


R m „(x) = 


ao + Six + 


+ x 

m 


m 


mn 


b 0 + bix + ••• + b n x 


(i.D 


is roughly equivalent in its "curve-fitting" ability to a polynomial of 
degree m+n. In some instances, however, such rational forms are far 
superior to polynomials. 

Heuristically , the consideration of rational functions as approxi- 
mating forms is motivated by a comparison with polynomial approximation. 
Consider the approximation of a function, f, of a single independent 
variable, x, by a polynomial: 


1 


2 


f(x) — Co + CiX + ••• + c n x n (1.2) 

The parameters to be determined, the c^, enter the problem in a linear 
fashion. For a rational form, however, the parameters, enter in a non- 
linear manner since 



a 0 + aix + ••• + a m x 1 " 
bo + bix + + b x n 


can be written as 

R m „(x) = r 0 + rix + ••• + r x m 
mn' ' m 

where 


(1.3) 


(1.4) 


r = (1.5) 

b 0 + bjx +<-■>» + b n x 

Unlike the c^, the r k are not constant but vary with x. Hence one would 
expect greater flexibility in approximating with a rational form than 
with a polynomial. In fact, a polynomial is but a special case of a 
rational form, for if 

b 0 + bix + ««• + b x n = 1 

then r^ = a k and one obtains the form of (1.4). 

Algorithms for. determining "best" rational approximations (optimum 
in some sense) have been found by various investigators (5), (6), (13). 
However, for optimality under the Tchebycheff norm these are, at best, 
involved and are usually iterative in nature. This investigation is 
concerned with the problem of approximating certain non-linear functions 
in a manner which allows greater facility in their handling. In particular, 
interest is in obtaining rational forms which approach optimality under 
the Tchebycheff norm in their approximation of such functions. 



3 


1.2 Background Study of Rational Approximation 

One of the earliest successful attempts to obtain a method of 

analytically developing rational approximations was due to H. Padd in 

1892. It is a simple but effective method based upon a series expansion 

of the function approximated. It suffers from the disadvantage of the 

» 

Taylor expansion in that for a finite ordered approximation, the error 
increases as one progresses further from the origin. In spite of this, 
the Pad6 method forms the basis of the method developed in this 
investigation. 

Apparently it was not until the late 1950's and early 1960‘s that 
rational approximation was extensively investigated. Shanks (17) in 
1954 investigated several useful classes of non-linear transformations 
which yield rational forms from both converging and diverging sequences. 
Shanks* efforts also proved that Padd approximation is but a -special 
case of his transformations. 

Wynn (21) examined the rational approximation of functions defined 
by a power series and developed the so-called "epsilon algorithm," also 
a special case of Shanks' non-linear transformations, 

Cheney (4) and Boehm (3) among others have developed considerable 
formal theory, having investigated existence, characterization, and 
convergence properties of rational Tchebycheff approximations (i.e., 
rational approximations which are optimal under the Tchebycheff norm). 

Both Cheney and Maehly (13) have developed a number of iterative algo- 
rithms for obtaining rational approximations which are optimal in the 
Tchebycheff sense. The work of these investigators is generally repre- 
sentative of the current level of development of optimal rational 
approximation. 
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1.3 Scope of the Investigation 

For this investigation a restriction is introduced regarding the 
nature of the function to be approximated. Each function, f, must be 
expandable in a Taylor series. 



The reason for this restriction rests with the Padd method which forms 
the basis for the development. 

An important distinction should be made. The concern here is 
with obtaining rational forms as approximations to functions expressed 
as power series. The investigation is not concerned with approximating 
a function knowing only a set of its values. 

Within the restrictions set forth above, the purpose of the 
investigation may be summarized. The purpose is to develop a relatively 
simple method for obtaining rational forms as approximations to functions 
admitting a power series representation, and which are near-optimal in 
the Tchebycheff sense in their approximations of such functions. The 
motivating question is, "To what extent can such functions be accurately 
approximated by rational forms?" 



CHAPTER 2 


APPROXIMATION CHARACTERISTICS 


This chapter sets forth some of the basic formal theory associated 
with the development of best rational approximations and is based on 
Cheney (4). 


2.1 Existence of Best Rational Approximations 

A family of rational functions, r , is first defined. 

mn 

r mn = { R mn E {j* = m, 3Q = n, Q(x) > 0 on [a,b]J- (2.1) 

r mn 1S c ^ ass national functions defined on the closed interval 
[a,b] where 


p = P m (x) = a 0 + aix + ••• + a x m 
m m 


( 2 . 2 ) 


Q = Q n (x) h b 0 + bix + + b n x n 

The polynomials P and Q are regular polynomials in the single independent 
variable, x. The degree (8) of P and Q satisfy the inequalities 

8P = m, 8Q = n (2.3) • 


where m and. n are integers specifying the order of the rational form, 
R mn‘ The i nec i ua li t ' ies (2.3) admit the possibility that the polynomials 
comprising R mn may have degrees less than m or n. The reason for this 
will be illustrated in the next chapter. P and Q are further restricted 
to have no factors in common other than constants. R = P/Q is then said 
to be in irreducible form. The members of r must be bounded. Thus it 
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is both necessary and sufficient that Q have no roots on the interval 
[a,b]. This may be accomplished without loss of generality by requiring 
that Q(x) > 0. 

The condition of optimality will be the minimization of the 
Tchebycheff norm (uniform norm, infinity norm), 

min || || = min f max | |l (2.4) 

lxe[a,b] > 

Defining the error function (f is the function approximated), 

S(x) s f(x) - R mn (x) (2.5) 

the scalar, X, is obtained as 

X = min || <5{x) || = min / max J S(x) | 1 ‘ J2.6) ■ 

\xC[a,b] ' 

Thus for R to be a best approximation the local minima and maxima of 
6(x) for all xC[a,b] must have the same absolute value, X (see Figure 
2 . 1 ). 



Figure 2.1 Condition of Optimality 
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Having defined the class r , and having set forth the condition 
of .optimality, an existence theorem is now given. This is due to Cheney 
(4) and is stated here without proof. 

Existence Theorem : To each function, f, continuous on the 

closed interval [a,b], there corresponds at least one best 
rational approximation from the class r mn [a,b]. 

2.2 Characterization of Best Rational Approximations 

The condition of minimum Tchebycheff norm is now extended to 
more fully characterize best rational approximations. The concept of 
T-alternations is introduced by first stating that there exist k values 
for x, Xi < x 2 < < x^, such that 6(x.) <$(x.. + ^) < 0, i = 1, 2, 

k-1. (In Figure 2.1 there are six values of x for which the condition 
holds.) Thus 6(x) changes sign k-1 times and is said to have k T- 
alternations. The condition of optimality may now be recast into the 
following criterion: 

Cri teri on : In order for the irreducible rational function 

R = P/Q,£ T mn to be a best approximation to the continuous 
function, f, it is necessary and sufficient that the error 
function, 6, have at least 2 + max j(m + 3Q) , (n + 3P)j> 
T-alternations. 

The reason for having to select the maximum of the bracketed 
terms is because of the inequalities (2.3). This is the criterion used 
in this investigation for judging best rational approximations. 

In addition to proving the above, Cheney also proves the 
following: 

Uniqueness Theorem : Best approximations in r mn [a,b] are always 

unique. 
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The representation of best rational approximations, however, is not unique 
since multiplication of the numerator and denominator polynomials by an 
arbitrary constant implies an infinite number of representations. 



CHAPTER 3 


PADE APPROXIMATION 


The Padg method of obtaining rational approximations to functions 
defined by a power series is presented in this chapter. A matrix formu- 
lation is utilized and the method is illustrated with a simple example. 

An expression for the error is also derived.. 


3.1 Development of the Pad<§ Method 

Consider a function, f, which may be expanded in a Taylor series 
about the origin. 

03 

f(x) = S c -j x 1 (3.1) 

i=0 1 


A rational function of the form (1.1) is desired which will 
approximate f(x), 

m 

Ea, x k 

1 = k „°- •- + S(x) (3.2) 

2 b. x J * 
j=0 3 

where 6(x) is the error in the approximation. Multiplying (3.2) by the 
denominator on the right gives 




i=0 




(3.3) 


9 
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which can be rewritten as 

bj bj (3.4) 

with = a^ = ••♦ = 0, and c^ = 0 for k < j. Neglecting the error 
term, the coefficients of the first n + m powers of x may be equated to 
zero to yield a system of n + m + 1 equations in n + m + 2 unknowns. 

This under- determi nance may be relieved by normalizing one of the b's, 
say b 0 > to 1. The result is then a system of n + m + 1 equations in 
n + m. + 1 unknowns. For clarity, assume for the illustration that m * n. 
The resulting equations are 

Co ” So 
Ci + Co bi = di 

c 2 ^ Ci bi + Co b 2 = ct 2 



c + ••• + ci b n . + c 0 b„ = a n 
n n-1 n n 


(3.5) 


Co + ••• + c b ,+c 
2n n+1 n-1 n 

The system (3.5) may be rearranged and 

notation: 


b n =0 

written in the following matrix 


1 0 ••• 0 

0 1 ••• 0 
* 

9 * • 

• \ • 
• » • 
« 

0 1 

0 ..... 0 

• • 
• • 
• • 

o ..... 0 



(3.6) 
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In general the size of the square coefficient matrix is (m+n+1) x 
(m+n+1); the vector of unknown a's and b's, and the vector of c‘s are 
both (tn+n+1) x 1. The matrix equation (3.6) may be partitioned as 


(3.7) 


‘i ; H° 


f* *■ 

a 


3* cS 

h 

t 

JL 





i 

I 

| 





.° « G . 


b 


,9. 


where the sub-matrices are given below. 

[i] = (n+1) x (n+1) identity matrix, 
[0] = n x (n+1) zero matrix. 


W - 


0 

-Co 


0 

0 


[G] = 


r c n-l"‘ 

oo « 


-Cl 


~ c 2n-l* ° ” c n 


, (n+1) x n 


, n x n 


(3.8) 


(a) = 


So 


(h) = 


Co 


(n+1) x 1 


(b) = 


bi 


(g) - 


'n+1 


'2n 


n x 1 


Of course, in the general case the dimensions of these are 
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[i] , (m+1) x (m+1) ; [0] , n x (m+1) 

[H] , (m+1) x n ; [G] , n x n 

(a) , (h), (m+1) x 1 ; (b) , (g) , n x 1 

Because of the presence of the zero sub-matrix, the b coefficients may be 
solved independent of the a's as 

(b) = KT 1 (g) (3.9) 

Hence only an n x n system of linear equations need be solved simultane- 
ously. The a coefficients may then be determined immediately from 

(a) = (h) - [H] (b) (3.10) 

3.2 Example - Tan^U) 

The Padd method may be effectively illustrated using the inverse 
tangent function, 

f(x) = tan -1 (x) = x - y- + | — j- + ••• (3.11) 

Let the order of the desired rational approximation be arbitrarily 
specified as m = n = 2. Note that the expansion for tan _1 (x) is not 
strictly of the form of (3.1). Placing it in this form may be accom- 
plished in either of two ways. One is by factoring and making the change 
of variable, y = x 2 , 

f(x) = xg(y) » x^l - ^ + J*- - + •••) (3.12) 

and then approximating g(y). The second approach is merely to include 
the missing even powers of x with zero coefficients. This latter 
approach may be designated as placing the series in fundamental form . 

It is particularly convenient to use this approach in Chapter 5, and 
provides additional insight in this example. Thus in fundamental form 
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f(x) = 0 + x + (0)x 2 - y- + (Ojx 1 * - + ••• 

From this the Pad6 equations are found to be 






Co 

- ao 

Co 

= 0 



Cl 

+ 

Cob i 

** ai 

Ci 

= 1 

c 2 


Cjbi 


c 0 b 2 

“ a 2 

> C2 

= 0 

C 3 


c 2 bi 


Cib 2 

= 0 

C3 

= -1/3 

Co 

+ 

c 3 bi 


c 2 b 2 

= 0 

C4 

= 0 


(recall b 0 is normalized to 1). In matrix form, 


”100 0 o' 


ao 


Co 

010 -Co 0 


ai 


Cl 

0 0 1 -Ci -c 0 


a 2 

= 

C 2 

0 0 0 -c 2 -Ci 


bi 


C 3 

rsr 

0 

0 

1 

0 

CO 

1 

0 




b 2 


Cif 


Corresponding to (3.8) the sub-matrices are 



1 

O 

O 

l 


"o’ 


11 

i — 1 

a 

0 

0 

, (h) - 

1 



_-l 

0 _ 


_ 0 _ 


11 

1 — 1 

CD 

UTJ 

~ 0 

— 

-1 

> (g) = 

P 1/3 1 


1/3 

0 



1 

O 


Utilizing (3.9) and (3.10) the solution vectors are 


(b) - EG] -1 (g) = 


0 

1/3 


0 

-1 

0 


(3.13) 


(3.14) 


(3.16) 


(3.17) 


(a) = (h) - [H] { b ) = 
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Hence the resulting rational approximation is 

Ri 2 (x) = * tan _1 x (3.18) 

1 + 1/3 x 2 

Note that R i2 from (3.18) belongs to r 22 and that 3P * 1 even though it 
was desired that m = 2. Thus here is an example for which strict 
inequality of thd first of (2.3) holds. 

Figure 3.1 illustrates the behavior of the error curves associated 
with Ri 2 and with the Taylor expansion used to obtain R i2 * Obviously for 

i 

the given number of terms the accuracy has been improved.- This is the 
advantage of the Padd method. 

While its advantage is improved accuracy, the Padd method has an 
important similarity to the Taylor series. In fact, R m0 , obtained trivi- 
ally from the Padd method is just the Taylor expansion to m + 1 terms. 

Thus Padd approximation is the rational analogue to the Taylor series. 

An interesting example from Leibnitz is the evaluation of the ■ 
expansion of 4 tan^x) at x = 1. Thus, 

4 tan“ 1 (x) = 4(1 - 1/3 + 1/5 - 1/7 +•••)=* (3.19) 

Consider now the rational approximation R 44 U) - tan - 1 (x). Solving the 
corresponding Padd equations formed from the expansion (3.13) gives 


a 0 

= 0 

bo 

= i 

ai 

= 1 

bi 

= 0 

a 2 

= 0 ' 

b 2 

= .85714286 

a 3 

= .52380952 

b3 

= 0 

a 4 

= 0 

b* 

= 8.57142857 x 10 " 2 


for which the resultant rational form is 
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aix + a 3 x 3 

R 34 (x) = * tan -1 (x) (3.19) 

1 + b 2 x 2 + b 4 x 4 

Again, 3P is less than the value sought, m = 4. Evaluating at x = 1, 
4 Rs4(1) = 3.1372549 

which yields an error of 4.3377 x 10” 3 . To obtain this accuracy using 

* 

the Taylor expansion requires 23, to 231 terms depending on summation 
schemes. Thus for many slowly convergent series, the Pad£ method 
provides a useful and accurate representation with little effort. 


3.3 Error Expression 

The basic equation from which the Padd coefficients are found is 
just the sum of the first (n+m+1) terms on the left of (3.4) set equal to 
zero: 



b j C k-j 



(3.20) 


Substituting this expression back into (3.4) and solving for the error, 
6, yields 


<5 (x) = 


oo n 

E E 

k=n+m+l j=0 


J k-j 



(3.21) 


For convergent series an approximation of the error may be obtained using 
the k = m+n+1 term in the numerator sum. 

(c . . i b 0 + c . bj + • •• + c b )x m+n+1 

<$( x )= \ m+n+1 -- m+n m+1 n / (3.22) 

bo + bjX + ••• + b^ x n 
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In the previous example for 4 R 3 «♦( 1) - 4 tan~ x (l) = 

4 ( Cg + c B bj + c 7 b 2 + c 6 b 3 + c 5 b 4 ) x 9 

6(1) 

1 + biX + b 2 x 2 + b 3 x 3 + b 4 x‘ t 

6(1) * J.lll - (.143) (.857) » .2(8.57 x 10" 2 ) 

L 1 + .857 + 8.57 x 10~ 2 

« .01028 
i 

This is admittedly a crude approximation to the actual error, 
few. more terms would’ certainly improve the value. 


(3.23) 


Taking a 



CHAPTER 4 


THE TAU METHOD 


This chapter reviews some properties of Tchebycheff polynomials, 
and then presents Lanczos' Tau method by a simple illustration employing 
a series solution to a simple linear differential equation - . 


4.1 Tchebycheff Polynomials 

All the common orthogonal polynomials may be derived as solutions 
to the linear, second order differential equation due to Gauss, 

x( 1 ~ x)y" + fy - (a + 8 + l)x]y' - aSy = 0 (4.1) 

The solution to this equation is the hypergeometri c function, 

F(a, 8, Y, x) - 1 + yry x + x 


(4.2) 


+ a(a+l)(af2)8(g+l)(8+2) . ... 
Y(y+lT( Y+2) *1*2*3 


where a, 8 and y are constants. The series is convergent for |x| < 1 
and is divergent for |x| > 1 except if the series terminates after a 
finite number of terms. When the constants are chosen as 
Y « 1/2 


8 = n ) 

n: 

a = -n ; 


real, non-negative, integer 


and x is transformed to the new variable.. 


x ^ 

x 2 
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the hypergeometric function yields the Tchebycheff polynomials. 

F(-n, n, 1/2; ^ ) = T n (€) (4.3) 

The roots of these polynomials occur at 
C k = cos [( 2k— l)ir/2n] , 

and their spacing is such that 

X = min j max |T (S)|[ * min ||T (S)|| • 1 (4.4) 

(5 C [-1,1] " I 

In other words, recalling section 2.2, each possesses n+1 T-alternations. 
This is the property which makes the Tchebycheff polynomials useful in 
approximation, particularly in the following sections. 

4.2 The Tau Method 

Introduced in 1938 by C. Lanczos (11), the Tau method utilizes 
the uniform norm property of the Tchebycheff polynomials to improve the 
solution of various linear systems. The method is most easily presented 
using a series solution to a simple differential equation. 

Consider the following first order, linear differential equation: 
y' - y = 0 , y(0) = l (4.5) 

Although the solution is simply y = e x , a power series solution is 
assumed for the illustration as 

CO 

y(x) = c k x k (4.6) 

k=0 K 

Substituting into (4.5) and equating powers of x yields a system of 
coefficient equations which leads to the recurrence relation, 

c - C n-1 _ Co 

n n n! 
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Then, in the presence of the boundary condition, y (6) = 1, (4.7) pro- 
vides the solution to the differential equation (4.5), 

Y 2 Y n 

y(x) = 1 + x f 2 ~ + ••• + pr (4.8) 

the limit of which, as n is obviously e x . Since a series solution 
is usually exact only in the limit, practically speaking it cannot satisfy 
the differential equation exactly. Substituting (4.8) into (4.5) pro- 
duces an error of (x n /n!). Thus there will always be an error which may 
be made arbitrarily small but never zero. 

Now Lanczos reasoned that the error incurred by truncating the 
series solution could be countered by introducing another variable, t , 
into the differential equation. Thus the original equation, which was 
solved approximately, is perturbed to result in one which may be solved 
exactly. At the same time, judiciously choosing the way in which the 
new variable would be introduced could improve the accuracy of the series 
solution. 

Suppose the new variable, t is multiplied by a Tchebycheff 

n > 

polynomial of order n and the result added, to the right side of the 
differential equation. In effect this is an approximation of the error, 
and the coefficients which result actually appear to improve the accuracy, 
of the series solution. Returning to the example, this operation 
modifies (4.5) to yield 

y* - y = T n T n (4.9) 

Now an additional requirement is that the range of x must be 
known so that the Tchebycheff polynomials may be appropriately scaled. 
Hence for x c [a,b], 



(* 44 ^) 
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(4.10) 

For the example let xC [0,l]. Then (4.-10) yields the so-called "shifted" 
Tchebycheff polynomials, T^*. Letting n = 2, 

y‘ - y = i 2 t 2 * (4.11) 

and the assumed solution becomes a quadratic: 

y(x) = c 0 + Cix + c 2 x 2 (4.12) 

Substituting this into (4.11) and using T 2 * in its polynomial form 
(T 2 * = 8x 2 - 8x + 2)[l] the desired system of coefficient equations is 
obtained as 

Cj - Cq = T 2 

2c 2 “ c i = ~8 t2 (4* 13) 

-c 2 = 8 t 2 

The solution of this system yields 

Co = 1, Ci = c 2 ~ 8/9, t 2 = -1/9 

for which 

y(x) = 1 + 8/9x + 8/9x 2 , x C [0,l] (4.14) 

Comparison of the Tau. method with the corresponding series solu- 
tion is shown in Figure 4.1, giving the corresponding error curves. 

Note that gaining in accuracy toward the right side of the interval 
requires sacrificing some of the accuracy toward -the left. In spite of 
this, the net result is a gain in accuracy, under the Tchebycheff norm, 
over the interval. 

The foregoing illustration of the Tau method utilized a simple 
first order, linear differential equation for which one T-variable was 




used. In the case of higher order equations additional t's may be 
required. For example, a second order equation containing derivatives 
of only second order would require two x-variables. Also, Lanczos 
pointed out that the method is applicable to any system of linear 
equations, and use is made of th.is fact in the next chapter. 



CHAPTER 5 


TAU AND PADE METHODS COMBINED 

As a system of linear algebraic equations, the Padd coefficient 
equations are well-suited for application of the Tau method. This 
chapter develops the procedure, again using matrix notation, and derives 
the associated error expression. The nature of the method is examined 
as applied to various transcendental functions. Approximate solutions 
and representations are obtained for the classic transcendental equation 
of Kepler. 


5.1 Development of the Tau-Padd Equations 

Following the development of the Padd method, a function, f, 
expandable in a Taylor series about the origin is to be approximated 
by a rational function of the form (1.1). 



(5.1) 


where, as before, 6 is the error. The function is assumed to be expres- 
sable in the fundamental form introduced in Chapter 3. As before, 
multiplying by the denominator term and rewriting yields equation (3.4). 
Neglecting the error term and taking the first N terms allows (3.4) to 
be written as 


24 



25 


b. c k _. x k - y, a k x k = 0 , c k _^ = 0, k < j 

m (5.2) 

N is an integer, greater than n+m, whose value will be determined later. 

Now, since N > n+m, equating coefficients of like powers of x results in 

an overdetermined system. According to the Tau method, the procedure 

would be to introduce a T-variable as the coefficient of a Tchebycheff 

polynomial. Here, however, greater flexibility is available because any 

number of the c. 's are already available. Thus not just one but any 

number, l, of T-variables may be introduced. If z is determined as 

Z = N-m-n, then the Z T-variables may be added to (5.2) resulting in 



j=C 


b j = 


22 .. \ T k M 


X + 4 J 'V 

k k=n+m+l k k 


(5.3) 


The value of N is thus seen to be just n+m+£. The Tchebycheff polynomials 
must, of course, be scaled to the appropriate interval of x. Writing the 
Tchebycheff polynomials in their powers of x gives 

T k (x) = 



where the coefficients s kr are found according to the relation 


s = {-l) r i , (k ~r )! 2 (k " 2r "^ 

s k,2r 1 L) k-r (k-2r)l 71 6 

s k,2r+l = °» r = °> 1 > **• k / 2 


(5.5) 


In this notation (5.3) becomes 







X 


r 


(5.6) 
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The double sum on the right may be reversed if the upper limit on r is 
set to N. Then (5.6) may be rewritten to yield the basic equation for 
the Tau-Padd combination, 


N n 


m 


2 S b j c k-j xk * a k xk + X) T '- s, '“ x 


r (5.7) 


k kr 


k=0 j=0 


k=0 


r=0 k=n+m+l 


with the understanding that c^ . = 0 for k < j and s^ = 0 for k < r. 
If a step function of the indeces is defined as 


u(k,j) = 


0, k < j 

1 , k > j 


(5.8) 


then the restrictions on c^ . and s^ r may be expressed qualitatively in 
the basic equation: 


N n m N N 

EE n c k-j xk ■ E*k xk +E E x k s kr u(k,r) 


N N 


J k-0 

k=0 j=0 k=0 r=0 k=n+m+l 


(5.9) 


The use of the step function, u, is most advantageous from a programming 
standpoint. 

Following the line of development of the Padd method, the coeffi- 
cients of like powers of x may be equated to yield the necessary system 
of linear coefficient equations. These may also be written in matrix 
notation. As before, let m=n for clarity, and normalize the b 0 coeffi- 
cient to one. The resulting system is then given in Figure 5.1. Note 
that the coefficient matrix is essentially the same as for the Padd 
method except that now it is augmented by the inclusion of the scaled 
Tchebycheff coefficients. Correspondingly, the, vector of unknowns is 
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augmented by the addition of the x's. Partitioning the system yields, 
as before, 


(5.10) 


I \ H 
1 


a 


h 

0 ! G 

* 

b 


_g_ 


with the appropriate sub-matrices defined as follows: 

[i] = identity matrix, (n+1) x (n+1) 

[o] = zero matrix, (N-n) x (n+1) 


[H] = 


‘ 20 + 1,0 


“2n+l,l 


‘2n+l,n 


’N,0 


S N,1 


’N,n 


0 

■c 0 


0 

0 


‘Vi*** 


(5.11) 


, (n+1) x (N-n) 


[G] = 


$ 2n+l,n+l *** s N,n+l 


‘N,N 


•. P • o • * 

L n 


-Ci 


-c 


N-l 


-c, 


N-n 


, (N-n) x (N-n) 


(h) = 


(a) = 


Co 

• 

• 

• 

, (n+1) x 1-; 

(g) = 

c n+l 

♦ 

• 

• 

_ c n_ 

* 


C N 

T 2n+1 

a 0 



• 

• 

• 

♦ 

♦ 

, (n+1) x 1; 

(b) = 

t N 

a 

> 


bi 

n 



♦ 

• 

_ b n _ 


(N-n) x 1 



In the general case the dimensions of the sub-matrices are 
[i], (m+l) x (m+l) ; [o], (N-nt-1) x (m+l) 

[H], (irH-1) x (N-m-1); H. (N-m-l.) x (N-m-1) (5.12) 

(a), (h), (m+l) x 1; (b), (g), (N-m-1) x 1 

The corresponding expressions for the solution of the augmented system 
are identical with those for the Padd coefficients, 

(b) = [G] _1 (a), (a) = (h) - [h] ( b ) (5.13) 

Since the application of the Tau method adds a t-variables to 
the Padd equations, it is convenient to alter the notation slightly to 
reflect the number of t-variables added. This is done by merely adding 
a third subscript to the symbol [^ . Hence a rational approximation of 
order m,n, found from the Tau-Padd equations using i -e-variables is 
indicated by R mnr 

5.2 Examples 

The following examples serve to illustrate the Tau-Padd combina- 
tion and to bring out some important aspects of the method. 

f(x) - e ♦ This first example already has its Taylor expansion 
in the fundamental form introduced in Chapter 3. 

f(x) = e x = 1 + x + 1/2! x 2 + 1/31 x 3 + ••• (5.14) 

For the desired rational form let m = n = 2. Let the number of T-variables, 
£, be six. Thus, N = £ + m + n = 10. The choice of i is rather arbitrary, 
being suggested primarily by experience. For rational .forms of order (2,2) 
rarely will more than six r's be required, and four will often suffice. 

For this example let x£ [0,l]; then the scaled Tchebycheff coefficients 
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correspond to the coefficients of the shifted Tchebycheff polynomials, 

_ t -,\r k (k-r)l „(k-2r-l) , c -, c \ 

s kr ” k-r f rT 2 * {5 - 15 ' 


r = 0, 1, 2, •••, k 


Under these specifications the components of the partitioned augmented 
coefficient matrix are 



0 0 
-Co 0 , 


S 52 *** “Cl “Co 


3x8 


S53 **• $ 1 0 > 3 ” C 2 “ C 1 



8x8 


0 ••• • S10 >10 “C9 -c 8 


( 5 . 16 ) 


(h) T = ( c 0 Ci c 2 ) 

(g) T 3 ( C3 C4 cs Ce C7 c 8 C9 C10 ) 


Solving the system according to (5.13) yields the set of values listed 
below. For comparison the corresponding Padd coefficients are also 
gi ven . 

Tau-Pad£ Pad£ 

3o ~ 1.0000031 ao = 1.0 

ai = .54164234 aj = .5 

a 2 = 


a 2 =• .10792084 


.0833333 
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t 5 = 3.49986928 x 10“ 6 

Tg = 4.36506101 x 10“ 7 

t 7 = 3. 09266682 x 10“ 8 

t 8 = 1.55708614 x 10" 9 

T 9 = 5.66904584 x 10" 11 

tio= 1.19460073 x 10“ 12 
bo = 1.0 bo 

bi = -.45821125 b* 

b 2 = 6.50542644 x 10“ 2 b 2 

x a 0 + aix + a 2 x 2 

e ~ R226 = b 0 + bix + b 2 x* ’ x 


1.0 

-.5 

.0833333 

0,1 


(5.17) 


Figure 5.2 plots the error functions associated with the Taylor series 
(6y), the Padd method (6^) , and the Tau-Padd combination (<5^). For the 
Taylor and Radd approximations five terms of the expansion (5.14) were 
used. Although eleven terms were used to determine the coefficients in 
(5.17), the results require no more computational effort for evaluation 
of the rational form than for the standard Padd form. The maximum of 
the Tau-Padd error, 6^, is three orders of magnitude less than for <5p 
and 6-j.. Under the criterion for an optimal approximation (Chapter 2), 
the Tau-Padd approximation is nearly optimum. The error curve has 
2 + max (m + 8Q, n + 8P) =6 alternations which approach T-alternations: 

min |j$ T || = 3.09 x 10" 6 , max j|S T j} = 6.68 x 10~ 6 


f(x) = ftn(l+x) . As a second illustration, consider the series 
expansion for &n(l+x) whose fundamental form is 

f(x) = An(l+x) = 0 + x - 1/2 x 2 + 1/3 x 3 - 1/4 x 4 + ••• (5.18) 
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As before an approximation of order 2,2 for x C [0,l] will be obtained 
using six x-variables. Solution of the system of augmented coefficient 
equations yields the following values which are compared with the 
standard Pad6 method: 




Tau-Padd 





Padd 

a® 

= 

-2.05651975 

X 

10“ 5 

ao 

- 

0.0 

ai 

= 

1.0009656 



a! 

s: 

1.0 



.62730344 



a2 

= 

.5 

T5 

= 

-2.47728605 

X 

10" 5 




t 6 


-5.59521770 

X 

10" 6 




T 7 

= 

-1.62740556 

X 

10" 6 




t 8 

= 

-2.70964554 

X 

10“' 7 




t 9 


-3.27816790 

X 

10" 8 




TlO 

= 

-1.63798206 

X 

10" 9 




bo 

= 

1.0 



b 0 

= 

1.0 

bi 

= 

1.1344666 



bi 

= 

1.0 

b 2 

= 

.21541081 



Pz 

ST 

.16666667 


Figure 5.3 shows the various error functions, 6 , 6 , 6 T , for the approxi- 
mations to £n(l+x)^ In this example the Tau-Pad£ error behaves similar 
to an optimal error curve up to x ^ 0.8. Thereafter, however, it grows 
rapidly, much as the Padd error. Also there is not as extensive an 
advantage with the Tau-Pad£ combination in this case as in the previous 
example, 6 and 6 being closer. The Taylor error appears to grow 

r 1 

extremely rapidly. This behavior is, of course, inherent in the nature 
of the defining series. Consideration of the convergence of the series 
offers a possible explanation. Because the convergence is quite slow, 
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six T-variables may be not be enough to obtain the desired near-optimal 
behavior. Also, near the upper limit of the interval the convergence 
becomes extremely slow-, divergence occurring at x = 1.0. 

f(x) ~ 1/(1 -fr 2x + x 2 . An example is now offered which is 
actually an irreducible rational function. Now one would expect that 
both the Padd and Tau-Padd methods should produce the function exactly. 

If the indicated division is carried out, the following series results: 

f(x) = 1 - 2x + 3x 2 - 4x 3 + ••• (5.14) 

For nK), n=2, the Padd method gives back the rational form of f(x). In 
fact, in all cases where f(x) is itself an irreducible rational function 
of the form (1.1), and m and n have the same values for the approximation 
as for the function, m and n, the Padd method simply returns the original 
function. Similarly, using the Tau-Padd combination with any number of 
T-variables gives the exact form of f(x), the T-variables being zero. In 
other words, if the function to be approximated is itself a member of 
r mn then the Tau-Padd method yields the function exactly regardless of 
the number of T-variables which may be introduced. 

A natural question is to inquire about the performance of the 
method when m and n do not happen to coincide with the actual order, m, 
n, of the rational function, f. First if either or both m and n are less, 
the resulting approximation merely will be a lower order approximation. 

Now if m > m, n = n or if m = m, n > n the Tau-Padd method will yield the 
correct a's and b's with 

a m+l = "• = a m = 0 for m > N, 

b~ + ^ = ••• = b = 0 for n > n. 
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up through x N be considered even though it does not appear explicitly 
in the series. The most consistent way to insure this is through re- 
duction to fundamental form. 

Proceeding with the example, let &=4 and again x £ [0,l].- Then 
the Tau-Padd equations yield the following parameter values: 


a 0 = 

8.06609950 x 10’ 8 

b 0 = 1.0 

ai = 

.9998780 

bi = .34753262 

a 2 = 

.34783041 

b 2 = 1.0408238 

a 3 = 

.70477263 

b 3 = .29418535 

at, = 

.19018504 

b 4 = .17348426 


t 9 = 1.06674826 x 10" 7 

T X o - 2.89900133 x 10’ 8 

Tn = 3.10652367 x 10’ 9 

Tia « 1.30341767 x 10“ 10 ' 

The error functions for both the standard Padd and Tau-Padd methods are 
shown in Figure 5.4. 

An evaluation of ir may again be obtained, 

4 IW(19 = 3.141114136 

with an error of 4.775 x 10“ 4 . To obtain the same accuracy by merely 
summing the expansion requires 2090 terms. 


5.3 Error Expression 


The error expression for the Tau-Pad£ combination is found the 
same way as for the standard Padd method. Substituting the basic 
equation (5.7) into the general expression (3.4) “and solving for S(x) 


PRECEDING PAGE BLANK NOT FILMED 
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k < r) 
( 5 . 21 ) 


An interesting result is the limit of (5.16) as N »: 


6 (x) = Lim 6(x) 
N + ~ 



• (5.22) 


For convergent series a rough approximation of the order of magnitude of 
the error may be obtained by this expression. In general the maximum 
value of each occurs at the interval limits and is ±1.0. When x is 
normalized such that |x| £ [o,l] , (5.22) reduces to 



6J1) = 3,969 6 q 7 10 6 = 6.539 x 10 -6 

which is slightly less than the actual maximum error, 6.68 x 10~ 6 , 
Caution must be used when using this method, however, since it loses 
validity when the approximations are not close to the optimal. 
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5.4 Applications to Kepler's Equation 

A fundamental equation of astrodynamics is the transcendental 
relation between time and angular position of an orbiting body. 

For elliptic motion: M = E - e sin E (5.24a) 

For hyperbolic motion: M - e sinh H - H (5,24b) 

M is the mean anomaly (a function of the mean motion and the time), e is 
the eccentricity, and E and H are the eccentric and hyperbolic anomalies, 
respectively. As transcendental equations, (5.24a,b) generally are 
solved for E and H using some iterative technique such as the method of 
successive approximations or the Newton-Raphson method. Such algorithms 
require reasonable starting values which are obtained by some method of 
approximation. In this section attention is given to the problem of 
obtaining solutions for E or H using the Tau-Pad6 equations. In the 
final paragraphs, Battin's modification (2) of Herrick’s universal var- 
iable formulation of Kepler's equation will be introduced, and a different 
application of the Tau-Pad£ combination will be made. 

Approximations to Sin and Sinh . • If rational approximations for 
the sine and hyperbolic sine functions are obtained in the form 

ao + aix + a 2 x 2 

R2i(x) mix — < 5 - 25 > 

then (5.24a) and (5.24b) can be written in the common form, 

( ao + aix + a 2 x 2 \ 

— ) (5.26) 

1 + bix / 

where x is E or H as required. Further, it is necessary to obtain the 
approximations only over the positive values of x, since both functions 
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are odd, and the sign on x is positive or negative according to the sign 
on M. The identification of elliptic or hyperbolic motion may be made 
depending on whether e < 1 or e > 1, respectively (the case e = 1 
represents parabolic motion; the time-position relation for this motion 
requires no iterative technique for solution). Thus (5.26) may be re- 
written as 


yM = x 


ere 


a 0 + ajx + a 2 x J 

1 + bjX 


(5.27) 


where 

a = M/|Mj , y = (l-e)/|l-e| 

Multiplying (5.27) by 1 + bjx and combining coefficients of like powers 
of x gives the quadratic, 

Ax 2 + Bx + C » 0 (5.28) 

where 

A = bi - oea 2 

B = 1 - aeai - yMb 

% 

C = -creao - yM 

which may then be solved for x. 

The values of the a and b coefficients will now be found for the 
trigonometric approximations. First, scaling of the Tchebycheff coef- 
ficients is done for x C [o,TrJ since sine is periodic and odd. For sinh 
the interval x C [o,tt/ 2] is chosen with the value tt/ 2 being selected 
arbitrarily since the hyperbolic. sine is not periodic. The number of 
x-variables selected is eight. The results are summarized below. 



42 


Sin X = x - |r + fr - yr + 

a 0 « -2.738634 x 10" 2 b 0 = 1.0 

a x = 1.271591 bi = -1.155438 x 10'“ 

a 2 = -0.404789 

sinh x = x + -St + + ••• 

O. 0. /• 


a 0 = 2.565447 x 10" 3 b 0 - 1.0 

a 3 = 0.948262 b 3 = -0.356916 

a 2 - -0.193801 

The plots of the error curves for these approximations are shown in 
Figures 5.5a,b. Note that the errors are relatively large. One reason 
is because the primary contributor to the approximations is the first 
term in the Taylor expansions. Another important effect is that gen- 
erally, the larger the interval, the harder it is to approximate the 
function. In view of this, one might expect that the resulting solutions 
of (5.28) would not yield very accurate values of x. This is indeed the 
case as seen in Figures 5.6a,b.* The question as to which root of the 
quadratic to use has not been formally answered. However, experience 
has shown that using 


x 

u 


(-B + u/ B 2 - 4AC ) /2A, 


! +l, ellipse 
-1, hyperbola 


Undoubtedly a betteA appAoxlmatlon Mould a esutt li one fiound 
R321 &oa &ln and &lnh, and Aeca& t [ 5 . 28 ) as a cubic. HoMeveA, cubic 
equations oa e not paAtlculaAly convenient to &olve. 
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Figure 5.6a Error Curves for E 



Figure 5.6b Error Curves for H 
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consistently yields compatible results. 

In Figures 5.6a,b there is another effect present which is due 
to the nature of Kepler's equation itself. As e •+■ 1 and x ->0, a region 
of non-uniform behavior is approached, emphasized by the fact that the 
error curves go off scale. Moulton (14) has investigated this for the 
elliptic case, and Russell (16) has shown, using perturbation methods, 
that the non-uniform region is approached when 1 - e cos E ->0. One 
empirical way of countering this is to set x - M when 
1 - e cos M < e ls elliptic motion, 
e cosh M - 1 < e 2 , hyperbolic motion, 
where ei and e 2 are some judiciously chosen values on the order of .01. 


Variable Transformation . A second way to improve the accuracy 
for elliptic motion is one motivated by a suggestion by Gottlieb (8). 
Define a new variable, 

w = E - M = e sin E (5.30) 

The limiting values are obviously -1 < w < 1. Kepler's equation may be 
rewritten in terms of w as 

w - e sin (w + M) (5.31) 


or 


w - e(sin w cos M + cos w sin M) (5.32) 

Now since the interval of interest is now smaller, one could expect better 
accuracy. However, to retain the quadratic form (5.28) polynomials must 
be used to approximate sin w and cos w which partially offsets the advan- 
tage of a smaller interval. Approximating sine and cosine by R 20 8 , 
sin w = So + SiW + s 2 w 2 
COS W = C 0 + CjW + c 2 w 2 


(5.33) 
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and defining 

Sjjj = sin M, = cos M (5.33) 

(5.32) becomes 

w = e C M (s 0 + SiVi + s 2 w 2 ) + e S M (c 0 + CiW + c 2 w 2 ) 
or 

(e C M s 2 + e S M c 2 )w 2 + ( e C M Si + e c a - l)w 

+ (e So + e S M Co) = 0 (5.35) 

which may then be solved for w, and E subsequently found. Recalling 
that sine is an odd function, and cosine, even, the shifted Tchebycheff 
polynomial coefficients may be' used in the Tau-Padd equations. Then the 
sine approximations must be multiplied by the o defined in (5.27) so 
that (5.35) becomes 

Aw 2 + B w + C = 0, (5.36) 

with 

A = e(C M a s 2 + S M c 2 ) 

^ — s o s i + e c i ~ 1 
C ~ e ( a So + c o ) 

Here, the negative sign is used with the radical in the expression for 
the root, i.e., u = -1 in (5.29). Finally, solution of the corresponding 
Tau-Padd equations for coefficients of the R 2 o8 approximations yields 

So = -4.63945315 x 10" 3 , c 0 = 1.0021690 

$a = 1.0851999, Ci = -3.48778236 x 10“ 2 

s 2 — - .23475756, c 2 = - .42972092 

The associated error curves are shown in Figure 5.7, exhibiting the near- 
optimal characteristics (maxima of approximately equal magnitude, and 
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four T-alternatlons). The reason that the sine error is larger than the 
cosine error Is due to the number of terms employed in the Taylor series, 
only one term being taken for the sine as opposed to two for the cosine. 
Now even though the sine and cosine functions are approximated by R 20 8 
instead of a higher order the advantage of a smaller interval of 

approximation predominates. This is seen in Figure 5.8, showing the 
improved error behavior for E with e = .5 and e = .99. Immediately 
obvious, however, is the irregular behavior of these error curves. This 
is due to the fact that w is formed as a combination of the two quadratic 
approximations for sine and cosine, as shown by (5.32). This in turn 
maps into the solution of (5.36) in a non-linear fashion and hence the 
optimality exhibited in Figure 5.7 is not present in Figure 5.8. 

Transcendental Functions for the Universal Variable Formulation . 
To this point, concern has been with obtaining approximate solutions of 
Kepler's equation for the eccentric or hyperbolic anomaly. A different 
application of the Tau-Padd combination is now made to a universal var- 
iable formulation of Kepler's equation. Such a form, originally intro- 
duced by S. Herrick, allows writing just one time — angular position 

t 

relationship for both elliptic and hyperbolic motion. A modified 
formulation is given by Battin (2) as 


fit = 


ro * Vo 

fi 


x 2 C(a 0 x 2 ) + (1 - r 0 ct Q )x 3 S(a 0 x 2 ) + r 0 x 

(5.37) 


where 
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7q, Vo = Initial position and velocity vectors, 
t - time from initial to final position, 
p = gravitational parameter of central body, 
a 0 ** 2/r 0 - v 0 2 /p J r 0 = |r 0 |, v 0 = \ v 0 1 

The parameter, x, is termed the universal variable, and defined as 

E - E 0 

j , elliptic motion 

x s | ^ (5.38) 

/ H - Ho 

\ , hyperbolic motion 


with E 0 , E, H 0 > H being the initial and final eccentric and hyperbolic 
anomalies respectively. The C and S are transcendental functions defined 
by Battin (2) as 


S(u) 



(5.39) 



(5.40) 


While these series converge rapidly, over large intervals the development 
of their Tau-Padd approximants to high order is useful and informative. 
First, let u = ct 0 x 2 C ^-(2ir) 2 , (2ir) 2 J. Here, no symmetry, of the S(u) 
and C(u) exists. Hence no advantage of symmetry is present, and the 
full range must be used. This presents a problem because the interval 
is Targe, and scaling of the Tchebycheff coefficients requires multiplying 
each coefficient by some power of the interval length. This quickly leads 
to numerical difficulties for high order approximations since, for example, 
the twelfth order Tchebycheff polynomial has a scaled coefficient of x 10 
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on the order of IQ 16 , while the lowest order coefficient is 1. This 
problem (which can be encountered often for large intervals) is readily 
circumvented by the alternative of scaling u to [-l,lj. Hence, 

u n 5 (40w) n (5.41) 


so that 



(5.42) 


C(w) 


1 40w , (40w) 2 

2T " ?T” + 5! 


(5.43) 


for w C With these changes, solution of the Tau-Pade equations 

for the coefficients of R 448 as approximations for C and S gives 


C: 


do = ,50000005 

a x = -1.3329425 . 

a 2 = 1.2170535 
a 3 = - .44140671 
a 4 = 5.73370383 x 10 
b 0 = 1.0 
b x = .66744544 

b 2 = .21448521 

b 3 = 4.03536694 x 10" 2 


t 9 = -1.67491428 x 10" 7 
T i0 = 4.70365685 x 10" 3 
Tn = -6.53462358 x 10" 9 
t i 2 = 5.91614229 x 10" 10 
Xu = -3.89852249 x 10" 11 
t i 4 = 1.98790997 x 10" 12 
tjs = -8.01643151 x 10" 14 
Tie = 2.72493162 x 10 _1S 


3.82170297 x 10 
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S: 


a 0 = .16666667 

ai = ~ .23582030 

a 2 = . 14947671 

a 3 = -4.03986018 x‘ 10" 2 

a if = 4.16704687 x Kf 3 

b 0 = 1.0 

bi = .58507792 

b 2 = .16225461 

b 3 = 2.58883507 x 10~ 2 

bit = 2.04894379 x 10“ 3 


T 9 

=: 

-5.63359592 

X 

10" 9 

TlO 


1.42027615 

X 

10' 9 

Tn 

= 

-1.78816313 

X 

10“ 10 

Tl 2 

s: 

1.47893824 

X 

10“ 11 

Tl3 

= 

-8.96491447 

X 

10" 13 

f 1 4 

= 

4.22974700 

X 

10“ 14 

t I5 

= 

-1.58967531 

X 

10" 1S 

Tie 

s 

5.04612583 

X 

10" 17 


The associated error curves for w C [~l,l]are shown in Figures 5.9a,b. 
Because of the rapid convergence of the series, the Padd and Taylor 
results are practically the same. Hence the two errors are shown as 


one. 








CHAPTER 6 


EXPLICIT FORMS AND GENERALIZATIONS 

The development of some explicit expressions for the coefficients 
and x-variables for Tau-Pad6 approximations is outlined in this chapter. 
The resulting expressions are, in part, recursive in nature thus pro- 
viding compact and easily handled formulas. Generalizations of the Pad£ 
and Tau-Padd algorithms are discussed briefly, noting that their utility 
is limited. 

6.1 Explicit Expressions 

Except for the simplist of problems, a computer must be used for 
solution of the Tau-Pad£ equations. With this in mind, efficiency may 
be gained in terms of computer time and storage by employing explicit 
formulas for some of the Tau-Padd parameters. The development of such 
formulas is illustrated by deriving the explicit forms for Rn^- The 
extension to higher order forms will be obvious; explicit formulas for 
the parameters of R^ and ^ are given in the Appendix. 


For Ru 4 the Tau-Pad£ equations are 

Co = 3 o + T3S30 + T4S40 t T5S50 + T 6 $ 6 0 ( 3 ) 

Ci + Cob 2 = ai + T3S31 + T£,s 41 + T5S51 + t 6 s 6 i (b) 

C2 t Cibi = T3S32 t T4S42 f T5S52 *F TgSg 2 (c) 

< 

C 3 + C2bi = T3S33 + T4S43 + T5S53 + TgSg 3 (d) ( 6 . 1 ) 

C4 + C3bi - T4S44 + T5S54 + TgSgij (e) 

C 5 + C4b x = T5S5 5 + TgSg 5 (f) 

Cg + Csb! = TgS 6 g (g) 
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As before, sj^ is the coefficient of x r in 1 ^. The approach for solving 
these equations here is to use Gauss reduction. Hence ( 6 . 1 g) is solved 
for x6 j 

xe = — - (c 6 + Csbi) (6.2a) 


Substituting into (6. If) x 5 is found as 


t 5 


.-U 

s 5 5 | 


S 6 5 S 5 5 

Cs - 7— c G + (c 4 - ~ ~ c 5) b 


■} 


S "' a \ ft 

66 S G6 

Similarly, expressions for the other x's are found. 


(6.2b) 


X 4 = 


1 


s 44 


S54 S 5 4 S $ 5 S $4 

1O4 “ C5 "f" ( ~ ~ ™ T ) C 5 


S 55 


S55S66 S 6 6 


t 3 


* [■■- 

-*{■ 


S54 S54SG5 S$4 

C4 + ^ ) C5 


S 55 


S55S66 S 6 6 


H 


S43 

S33 f C3 ~ 


C4 + ( 


S43S54 S53 


S44S55 S55 


)cs 


S43S54S65 S43S64 S53S65 S 5 3 

+ ( - -I + + 


[ 


S44S55S66 

S 43 


S44S66 S55S66 $66 


) c 6 


c 2 - 


S 44 


C 3 + ( 


S 4 3S 54S 6 5 S53 


S44S55S66 S 55 


.) c 4 


S 4 3 S S 4 S 6 5 S43S64 S53S65 

+ ( _ — _ + + 


S 6 3 


S44S5SS6G S44S66 S55S66 Ss 6 


) C.] bj j 


(6.2c) 


(6. 2d) 


Substituting these expressions into ( 6 . 1 c) yields 
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S 32 S 42 S 32 S 43 

-C2 + - — C 3 + ( r 7 — 7 — ) c 4 + ( 


S 32 S 4 3 S 54 S 32 S 5 3 S 42 S 54 


$33 ’ ' S 44 S 33 S 44 * "* ' S 33 S 44 S 55 S 3 3 S 5 5 S 4 4 S 5 5 


S 52 . , S 32 S 43 S 54 S 65 S 32 S 43 S 64 . S 32 S 53 S 65 S 32 S 63 

*j- - ) C *f f «■ — — — ■ — — — — . — — — 4' - 1 ■ — ■ — - — ■ ■ — — * — < — — — — 

$55 $3 3 S 4 4 S 5 5 S 6 6 S 33 S 44 S 66 S 3 3 S 5 5 S 6 6 S 3 3 S 6 6 


S 42 S 5 4 S 6 5 


$ G u$ 


64^42 ^52^65 


S C oS I 


$44 S 55 S 66 


$ R G$ 


66^44 


$5 5$ 


5 5^66 


s 6 2 . r 

— ) C 6 + J- C l 

b 6S L 


'32 


+ C 

33 


S - - ~ 2 


(6.2e) 


+ (; 


'42 


*44 


S 32 S 43 
S 3 3 S 4 4 


) C 3 + ( 


S 32 S 4 3 S 54 


S 3 2 S 5 3 


S 42 S 54 S 52 


S 33 S 44 S 55 S 33 S 55 


S 44 S 55 S 55 


) C t 


S 32 S 43 S 54 S 65 S 32 S 43 S 64 S 32 S 53 S 65 S 32 S 63 S 42 S 54S65 

+ ( + — + ; + 


■ T — — l ’ - - — i 

S 33 S 44Ss5 S 66 $33 S 44 S 66 S 33 S5 5 S 66 S 33 Sgg S 44 S 55 S l 


S 64^42 

S 66 S 44 


S 52 S 65 S 62 

s 55 s 66 s ee 



0 


This expression may then be solved for b 2 . However, in its present 
form (6.2e) is extremely cumbersome. Examination of the terms yields 
a convenience which greatly facilitates the determination of the unknown 
coefficients. Judiciously grouping the rational forms of the s's 
results in the following recursion formulas: 

S 32 


1 

D 2 “ 7 ($42 ~ D 1 S 43 ) 

S44 

( 6 . 3 ) 

D3 = (S 52 “ 0iS$3 - D2S54) 

S55 



($6 2 " Dl s 63 " D2 s 64 ~ D 3 S 6S ) 
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Then in terms of these quantities, b 3 is simply 


~c 2 + Dic 3 + D2C4 + D3C5 + D 4 c 6 

-Ci + DiC 2 + D2C3 + D3C4 + D4C5 

The x-variables are then easily obtained by using equations (6. Id) - 
(6.1g). Similarly, a 0 and ai may be found immediately from (6.1a) and 
(6.1b). As the. foregoing implies, the algebra involved in the develop- 
ments is quite tedious for even the simple forms of R mn ^. By induction 
on the formulas for Rio 4 » Rii 4 S R 214 J and R 2 2 4 » those for the general 



(6.4) 


case, R„ 


, were found and are given in the Appendix. 


6.2 Generalization 

In Chapter 3 the classical Padd method was presented. It was 
shown to be applicable to those functions possessing a Taylor expansion. 
In fact, the Pad£ method is applicable to any function which can be 
represented as a linear combination of certain functions, for example, 
Legendre polynomials. Cheney has shown this in (4). By way of 
illustration, let f(x) be represented as 


f(x) = * k (x) (6.5) 

k=0 


where the <j>, are functions of the single variable, x. A rational 
approximation to f. 


a 0 <j> 0 + ai <|>2 + ••• + a <j> 

R (6.6) 

b 0 <j> 0 + bj <j>! + ••• + b n 4> n 

may be obtained in a manner analogous to the standard Pade method as long 
as the ^ satisfy the relation 
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i+j 

N v2>i 

k=0 


jk fyc 


(6.7) 


where the are constant coefficients. Obviously this fi-s the Padd 
method in general form. The classical Pad£ method is but a special case 
of this since then the are just the x . Equation (6,7) is satisfied 
since 


x 


i 



K-k - k " i+ j 

x i+j =» 

( A ijk = x " k = 1+ J 


( 6 . 8 ) 


An interesting development concerning the generalized Padg methbd 
was made by Maehly and is briefly presented by Snyder (18). Here the 
<$>|, are Tchebycheff polynomials. The implementation of the algorithm is 
not a difficult task since the corresponding .form of (6.7) is relatively 
simple: 


■ 



(6.9)' 


Generalizing the Tau-Padd equations is similar to the generalized 
Pad6. However, in addition to the requirement that the <f^ satisfy (6.7), 
the Tchbycheff polynomials also must be expressed as a linear combination 
of the 

n 

T n<*> = Z ?nk *k < 6 ' 10 > 

k=0 


It is particularly interesting to consider the problem of obtaining the 
Tau-Padd approximant when f(x) is of the form (6.5) with <f> k = T^. In' 
this case the procedure is nothing more than the application of Maehly 's 
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Tchebycheff form of the generalized Padd method. The introduction of the 
T-variables has no advantageous effect on determining the coefficients 
for the rational form. Thus the implication is that the Tau-Pade\ method, 
in effect, attempts to convert a power series into a "quasi -Tchebycheff" 
series (it will never actually make it since an infinity of terms would 
be required) while at the same time making use of the flexibility of 
rational forms. 

As one might suspect, relations (6.7) and (6.10) are generally 
not. readily available for tf^’s other than the simple powers of x or the 
Tchebycheff polynomials. Even when they are, the algebraic manipulations 
become overwhelming as shown by the form of (6.7) when the ^ happen to 
be P^, the Legendre polynomials (20): 


p p _ ^m-r ^r \-r / 2n + 2m - 4r + 1 \ 
n ^ n \2n + 2m - 2r + 1 / 


m 


n+m-2r * 


r=0 


n+m-r 


( 6 . 11 ) 


A = 
m 


(2m - 1) 


m: 


For these reasons, application of the generalized Pade or Tau-Pade method 
is usually not particularly attractive. 

A natural question is to ask if the algorithm can be extended to 
handle functions of several variables. For example, consider a -function 
of two variables represented by the following series: 


co 

f (x,y) = ^ 


k~0 


( 6 . 12 ) 
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Since the c^'s are no longer constants, direct application .of the Tau- 
Padd algorithm results in a rational approximation of the form 


m m 


^r— ■ \ \ 

2^ L, ’i^ 1 ¥ x) 


R (v V ) = t=0 J =Q 

K mr\V x,y} n n 


2 E »,« .,w 

i-0 j=0 


( 6 . 13 ) 


This approach was attempted on the series associated with the gravita- 
tional potential of the Earth (zonal harmonics only), the corresponding 
form of (6.12) being 



where y is the gravitational parameter of the Earth, r is the distance - 
from the coordinate system origin, and <|) is the geocentric latitude. 
Unfortunately, the amount of effort involved is extensive and no mean- 
ingful results have been obtained with this approach. Further investi- 
gation is necessary to determine its utility. 

Probably one of the more notable generalizations of the Tau-Padd 
algorithm is the realization that a different error behavior may be 
obtained by merely using a set of weighting functions other than the 
Tchebycheff polynomials. Thus, using the Legendre polynomials, for 
example, would yield an approximation with different error characteristics. 



CHAPTER 7 


SUMMARY 

This thesis is concerned with the problem of obtaining approxi- 
mations to functions defined by a power series in a manner which allows 
facility in their handling and which yields near-optimality under the 
Tchebycheff norm. Rational functions are chosen as the approximating 
form because of their simplicity and flexibility, and are shown to be 
superior in many cases to the more common polynomial approximation. The- 
results of this investigation offer-.a method which yields just such 
approximations. 

The classical Padd method forms the basis of the technique 
developed in the investigation, utilizing the Taylor series representation 
of the function to be approximated. One of the earliest analytical 
schemes, it is also one of the simplest, merely requiring the solution 
of a system of linear equations. 

The Tau method developed by Lanczos provides a tool for modifying 
the Padd method. It employs the important uniform error properties of 
the Tchebycheff polynomials to weight the coefficients in the Padd 
rational form to obtain near-optimal behavior of the associated error 
function. The combination of these two methods forms what, in this 
thesis, is referred to as the Tau-Padd method. It is illustrated by 
the approximation of numerous transcendental functions, and is shown 
to return the approximated function exactly when that function is of 
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rational form. Applications to the classical equation of Kepler offer 
additional insight into the method. 

Some important conclusions are drawn concerning the use of the 
Tau-Padd method. The approach is found to be effective in obtaining a 
rational form which is near-optimal in its approximation of functions 
admitting a power series representation. Since polynomials are but 
special cases of rational forms, the method may be easily used to obtain 
near-optimal polynomial approximations. In generalizing the Tau-Padd 
method to functions defined by series other than power series, difficul- 
ties are usually encountered which hamper implementation of the method. 
However, when functions are represented by a Tchebycheff series, the 
generalized Tau-Padd equations reduce -to a generalized Padd algorithm. 
From this the conclusion is drawn that the regular Tau-Padd method, in 
effect, attempts to "convert" a defining power series into a -Tchebycheff 
series. 

Two important restrictions are placed upon the Tau-Padd approach 
First, the function to be approximated must, of course, be expressible 
in the form of a known power series. Second, because of this, the algo- 
rithm obviously may not be employed to obtain rational approximations 
from discrete data values. The specification of an. interval over which 
the approximation is valid is an additional restriction but certainly 
not a severe one. 

In conclusion, the Tau-Padd method offers the capability of 
approximating functions defined by a power series. The technique pro- 
duces a near-optimal approximation in a non-iterative manner. It 
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provides an answer to the motivating question, "To what extent can such 
functions be accurately approximated?" 
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APPENDIX 


EXPLICIT FORMULAS FOR THE TAU-PADE COEFFICIENTS 



ao + ajx + • • • + a x m 
• 

bo + biX + • • • + b^ x n 


R mo&* 




k = m 


a 


k 



s m+j,k 


k < m 



j s 1, 4 ; (no sum for i > j-1) 



r = m+1, N; (no sum for i > N) 
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R mn£* 


n > 0: 


a,. = c, 


-r 

j=i 


T m+n+j s m+n+j ,k 


k = 0, 


m 


T 


r 






i=r+l 



r = m+n+1, N; (no sum for i > N) 



1 

s m+n+j ,m+n+ j 



j = 1, •••, l ; i = m+1, m+n; (no sum for k > j-l) 




c m+n+j-k 




k = 0 , • • • , n ; i = m+1 , • » » , m+n 


r 


b >l l W m+l,l Vl.2 * * * W m+!,n 


-i -1 


W. 


m+1,0 


W iW _ * • • W W ^ 

m+n,l m+n, 2 m+n,n m+n,0 


(n x 1) 


(n x n) 


(n x 1) 



